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' Full relativistic simulations in three dimensions are known to develop runaway modes 

that grow exponentially and are accompanied by violations of the Hamiltonian and 
momentum constraints. We present here a method that controls the violation of these 
Q ' constraints and is tested with simulations of binary neutron stars in circular orbits. We 

' show that this technique improves the overall quality of the simulations. 
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1. Introduction 



For most of the past decade, the main theoretical thrust in gravitational research has 
been directed toward obtaining stable and accurate numerical models of compact- 
object binary systems. One of the most difficult problems to tackle has been the 
control of exponentially growing instabilities that degrade the quality of any sim- 
ulation and, eventually, terminate it. We propose in Refs. 1,2 an evolution scheme 
, (Constraint relaxation or CR) where the Hamiltonian and momentum constraints 

are only approximately solved at every time step, gently steering the evolution to- 
ward the space of solutions of the Einstein Field Equations without completely 
forcing their numerical satisfaction. This method utilizes the conformal decompo- 
sition of the spatial metric and extrinsic curvature, which has traditionally been 
used to solve the initial value problem for binary systems. In this decomposition, a 
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conformal factor ip factored out of the spatial metric and a longitudinal addition to 
the extrinsic curvature generated from a vector potential Wi are used to satisfy the 
Hamiltonian and the three components of the momentum constraint respectively. 
rS ' CR drives tp to the solution space of the Hamiltonian constraint by means of a 

parabolic equation for the conformal factor. Similarly, the momentum constraint is 
controlled by the use of Wi to push the simulation toward the space of solutions of 
the momentum constraint. In both cases, a full relaxation of ip and Wi would lead 
to the numerical solution of the constraints. However, the stability of the relaxation 
methods relies on gently updating these fields during the evolution. We showed in 
Ref. 2 that a full relaxation scheme becomes unstable rather quickly when used in 
combination with BSSN. 

2. Results 

We tested the CR scheme with simulations of an irrotational BNS system. The 
details of the initial data set are provided in table 1 of Ref. 1 and the corresponding 
convergence tests are presented in Ref. 2 The plots show curves that, for clarity, 
have been normalized to their corresponding initial values. 

The evolution of the Li norm of the Hamiltonian constraint residual TLC and the 
three components of the momentum constraint residual Mi across the numerical 
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Fig. 1. Evolution of the L2 norm of the Hamiltonian constraint violation. 
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Fig. 2. Evolution of the L2 norm of the momentum constraint violation. 
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grid are shown in Figs. [T] and [2] respectively. CR achieves an impressive reduction 
of the Hamiltonian contraint violation (FigHJ). Note that the curves are plotted 
in logarithmic scale to highlight the more than two orders of magnitude difference 
between the CR (dashed) and BSSN (solid) results. CR not only suppresses the 
constraint violation modes, but also reduces the violation present in the ID set by 
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a factor of about 5. The difference between the CR (solid) and the BSSN (dashed) 
results in the case of the momentum constraint (Fig. [5]) is not that impressive. At 
the end of the simulation, the momentum constraint violation was about four times 
smaller than in the BSSN runs. The spikes present in the BSSN curves at t ~ 0.4P 
occur on the stellar surface and are related to matter displacement in the grid, a 
side-effect of using a frozen shift vector. Note, however, that those spikes disappear 
when using CR. 
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Figure [3] shows the evolution of the total angular momentum. The plot shows the 
CR (solid), BSSN (dashed) curves, as well as the PN estimation (dotted line) of the 
angular momentum loss for a point-mass binary with the same mass and angular 
momentum as the BNS in consideration (see Appendix B of Rcf. 1 ). The CR based 
run agrees with the PN prediction for about 1.5 orbital periods. The inset of figure 
[3] zooms in on the first half of the period, showing the reduced level of noise in the 
CR curve. 
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